lib/Mi manual de algoritmos/version_maraton_nacional_2008/src/geometria/is_inside_concave_polygon.cpp

   1 /*********
   2  * Point *
   3  *********
   4  * A simple point class used by some of the routines below.
   5  * Anything else that supports .x and .y will work just as
   6  * well. There are 2 variants - double and int.
   7  **/
   8 struct P { double x, y; P() {}; P( double q, double w ) : x( q ), y( w ) {} };
   9 struct P { int x, y; P() {}; P( int q, int w ) : x( q ), y( w ) {} };
  10
  11 /***************
  12  * Polar angle *
  13  ***************
  14  * Returns an angle in the range [0, 2*Pi) of a given Cartesian point.
  15  * If the point is (0,0), -1.0 is returned.
  16  * REQUIRES:
  17  * include math.h
  18  * define EPS 0.000000001, or your choice
  19  * P has members x and y.
  20  **/
  21 double polarAngle( P p )
  22 {
  23     if( fabs( p.x ) <= EPS && fabs( p.y ) <= EPS ) return -1.0;
  24     if( fabs( p.x ) <= EPS ) return ( p.y > EPS ? 1.0 : 3.0 ) * acos( 0 );
  25     double theta = atan( 1.0 * p.y / p.x );
  26     if( p.x > EPS ) return( p.y >= -EPS ? theta : ( 4 * acos( 0 ) + theta ) );
  27     return( 2 * acos( 0 ) + theta );
  28 }
  29
  30 /************************
  31  * Point inside polygon *
  32  ************************
  33  * Returns true iff p is inside poly.
  34  * PRE: The vertices of poly are ordered (either clockwise or
  35  *      counter-clockwise.
  36  * POST: Modify code inside to handle the special case of "on an edge".
  37  * REQUIRES:
  38  * polarAngle()
  39  * include math.h
  40  * include vector
  41  * define EPS 0.000000001, or your choice
  42  **/
  43 bool pointInPoly( P p, vector< P > &poly )
  44 {
  45     int n = poly.size();
  46     double ang = 0.0;
  47     for( int i = n - 1, j = 0; j < n; i = j++ )
  48     {
  49         P v( poly[i].x - p.x, poly[i].y - p.y );
  50         P w( poly[j].x - p.x, poly[j].y - p.y );
  51         double va = polarAngle( v );
  52         double wa = polarAngle( w );
  53         double xx = wa - va;
  54         if( va < -0.5 || wa < -0.5 || fabs( fabs( xx ) - 2 * acos( 0 ) ) < EPS )
  55         {
  56             // POINT IS ON THE EDGE
  57             assert( false );
  58             ang += 2 * acos( 0 );
  59             continue;
  60         }
  61         if( xx < -2 * acos( 0 ) ) ang += xx + 4 * acos( 0 );
  62         else if( xx > 2 * acos( 0 ) ) ang += xx - 4 * acos( 0 );
  63         else ang += xx;
  64     }
  65     return( ang * ang > 1.0 );
  66 }